Load Data

library(ggplot2)
library(data.table)
library(scales)
library(RColorBrewer)

data_dir <- here::here('..','data')
load(file.path(data_dir, 'exh.sampled.Rdata'))

car_set <- c('41BB','CD28','BAFF-R','CD40','TACI')


# map day 0 to both plus and minus
map_day_0 <- function(df) {
  return(rbind(
    df[k562!='none'],
    df[k562=='none'][, k562 := 'cd19+'],
    df[k562=='none'][, k562 := 'cd19-']
  ))
}

Assess Marker Thresholds

# table of thresholds chosen
marker_thresholds <- rbind(
    c(exh_opt_cd4$marker_thresholds, t_type='cd4'),
    c(exh_opt_cd8$marker_thresholds, t_type='cd8'))
marker_thresholds <- melt(
    data.table(marker_thresholds), id.vars=c('t_type'))
marker_thresholds[, value := unlist(value)]
marker_thresholds[, t_type := unlist(t_type)]
marker_thresholds <- rbind(
    copy(marker_thresholds[, k562 := 'cd19+']), 
    copy(marker_thresholds[, k562 := 'cd19-']))

blue_shades <- brewer.pal(7, "YlGnBu")[3:7]

# assess marker thresholds
ggplot(exh_dt[k562 %in% c('cd19+','cd19-') & variable != 'zombie_t' & 
    variable %in% names(exh_opt_cd4$marker_thresholds)]) + 
    geom_density(aes(x=value, color=factor(day))) + 
    scale_x_continuous(limits=c(-1000,5000)) + 
    theme_minimal() +
    geom_vline(data=marker_thresholds[variable != 'zombie_t'], aes(xintercept=value)) +
    facet_grid(t_type+k562~variable, scale='free') + 
    scale_color_manual(values=blue_shades)
## Warning: Removed 10287 rows containing non-finite values (stat_density).

Summary Stats by MFI, % Positive, and % Positive MFI

# which markers for which summary stats
total_mfi <- c('cd39_t','myc_t', 'pd1_t')
pos_mfi <- c('tim3_t','lag3_t','cd39_t','pd1_t')
pos_pct <- c('tim3_t','lag3_t','cd39_t', 'pd1_t')

pct_pos_markers <- exh_dt[variable %in% pos_pct,
  list(pct_pos= sum(gt_thresh)/.N),
  by=c('plate','well','day','k562','donor','car','variable','t_type')][, 
    list(mean_pct_pos= mean(pct_pos), sd_pct_pos= sd(pct_pos)),
    by=c('day','k562','donor','car','variable','t_type')]
pct_pos_markers <- map_day_0(pct_pos_markers)

pos_mfi_markers <- exh_dt[variable %in% pos_mfi,
  list(pos_mfi_well= mean(value[gt_thresh==T]), pos_mfi_well_sd= sd(value[gt_thresh==T])), 
  by=c('plate','well','day','k562','donor','car','variable','t_type')][,
    `:=`(mean_mfi_pos= mean(pos_mfi_well), sd_mfi_pos= sd(pos_mfi_well)),
        by=c('day','k562','donor','car','variable','t_type')]  
    
pct_pos_markers[, mean_diff_pct_pos := mean_pct_pos-mean(mean_pct_pos, na.rm=T), 
    by=c('day','k562','donor','variable','t_type')]

pct_pos_lbls <- pct_pos_markers[!is.na(variable) & !is.na(donor) & k562 == 'cd19+',
  list(var_mean = mean(mean_pct_pos, na.rm=T)), by=c('k562','day','variable','donor','t_type')][,
    variable := as.character(variable)]

# all CARs, CD19+ only

ggplot(pct_pos_markers[!is.na(variable) & !is.na(donor) & k562 == 'cd19+' & variable %in% pos_pct]) + 
    geom_bar(aes(x=car, y=mean_diff_pct_pos, fill=variable), stat='identity') +
    geom_linerange(aes(x=car, ymin=mean_diff_pct_pos-sd_pct_pos, ymax=mean_diff_pct_pos+sd_pct_pos)) +
    geom_text(data=pct_pos_lbls, 
        aes(label=percent(round(var_mean,3))), y=Inf, x=-Inf, hjust=0, vjust=1) + 
    facet_grid(variable+donor~day+t_type, scales='free') +
    theme_minimal(base_size=18) + scale_y_continuous(labels = percent_format(accuracy = 1)) +
    theme(axis.text.x=element_text(angle = 90, hjust = 1)) +
    labs(x='CAR', y='Change in % Positive', title='Pct % in Exhaustion Marker from Average CAR (All CARs)') +
    geom_hline(yintercept=0) +
    geom_vline(xintercept=-Inf)
## Warning: Removed 52 rows containing missing values (position_stack).
## Warning: Removed 52 rows containing missing values (geom_linerange).

# selected car set only, CD19+ only

pct_pos_markers[car %in% car_set, mean_diff_pct_pos := mean_pct_pos-mean(mean_pct_pos, na.rm=T), 
    by=c('day','k562','donor','variable','t_type')]

pct_pos_lbls <- pct_pos_markers[car %in% car_set & !is.na(variable) & !is.na(donor) & k562 == 'cd19+',
  list(var_mean = mean(mean_pct_pos, na.rm=T)), by=c('k562','day','variable','donor','t_type')][,
    variable := as.character(variable)]

ggplot(pct_pos_markers[car %in% car_set & !is.na(variable) & !is.na(donor) & k562 == 'cd19+' & variable %in% pos_pct]) + 
    geom_bar(aes(x=car, y=mean_diff_pct_pos, fill=variable), stat='identity') +
    geom_linerange(aes(x=car, ymin=mean_diff_pct_pos-sd_pct_pos, ymax=mean_diff_pct_pos+sd_pct_pos)) +
    geom_text(data=pct_pos_lbls, 
        aes(label=percent(round(var_mean,3))), y=Inf, x=-Inf, hjust=0, vjust=1) + 
    facet_grid(variable+donor~t_type + day, scales='free') +
    theme_minimal(base_size=18) + scale_y_continuous(labels = percent_format(accuracy = 1)) +
    theme(axis.text.x=element_text(angle = 90, hjust = 1)) +
    labs(x='CAR', y='Change in % Positive', title='Pct % in Exhaustion Marker from Average CAR (Selected CARs)') +
    geom_hline(yintercept=0) +
    geom_vline(xintercept=-Inf)
## Warning: Removed 20 rows containing missing values (position_stack).
## Warning: Removed 20 rows containing missing values (geom_linerange).

Number of positive markers per cell

marker_ct <- exh_dt[variable %in% c('tim3_t','lag3_t','cd39_t', 'pd1_t'),
  list(n_pos_marker= sum(gt_thresh)), 
  by=c('well','plate','day','Time','car','donor','k562','t_type')]

marker_ct <- marker_ct[, N := .N, by=c('well','plate','day','car','k562','t_type')][N > 1]

marker_ct <- map_day_0(marker_ct)

marker_freq <- marker_ct[, data.frame(table(factor(n_pos_marker, levels= 0:4))), 
  by=c('well','plate','day','car','donor','k562','t_type')]
names(marker_freq)[names(marker_freq) == 'Var1'] <- 'n_pos_marker'

marker_freq[ n_pos_marker != 7, pct := Freq/sum(Freq), 
  by=c('well','plate','day','car','donor','k562','t_type')]

marker_freq <- marker_freq[,
  list(mean_pct= mean(pct, na.rm=T)),
  by=c('day','car','donor','k562', 'n_pos_marker','t_type')]

ggplot(marker_freq[k562 != 'none' & n_pos_marker != 7]) + 
  geom_bar(aes(x=car, y=mean_pct, fill=n_pos_marker), stat='identity') + 
  facet_grid(k562 + donor + t_type ~ day, scales='free') + 
  scale_y_continuous(expand=c(0,0)) + 
  scale_fill_brewer('Number of Exh.\nMarkers Present', palette='YlOrRd') + 
  theme_minimal() +
  labs(x="CAR", y='Fraction of Cells')

ggplot(marker_freq[k562 != 'none' & n_pos_marker != 7 & k562 == 'cd19+']) + 
  geom_bar(aes(x=factor(day), y=mean_pct, fill=n_pos_marker), stat='identity', width=1) + 
  facet_grid(donor ~ t_type + car, scales='free') + 
  scale_y_continuous(expand=c(0,0), labels = percent_format(accuracy = 1)) + 
  scale_x_discrete(expand=c(0,0)) +
  scale_fill_brewer('Number of Exh.\nMarkers Present', palette='YlOrRd') + 
  theme_minimal() +
  geom_vline(xintercept=0:7-0.5) +
  geom_hline(yintercept=c(0,1)) +
  labs(x="CAR", y='Fraction of Cells')

Pie Charts, inidividual donors

ggplot(marker_freq[donor == 1 & k562 != 'none' & n_pos_marker != 7 & k562 == 'cd19+']) + 
    geom_bar(aes(x=1, y=mean_pct, fill=n_pos_marker), stat='identity', width=1) + 
    facet_grid(t_type + day ~ car) + 
    scale_y_continuous(expand=c(0,0)) + 
    scale_x_discrete(expand=c(0,0)) +
    scale_fill_brewer('Number of Exh.\nMarkers Present', palette='YlOrRd') + 
    theme_minimal() +
    theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
    geom_vline(xintercept=0:7-0.5) +
        labs(x="CD4/CD8 and Day", y='CAR', title='Fraction of Cells (Donor 1)') +
    coord_polar("y", start=0)

ggplot(marker_freq[donor == 2 & k562 != 'none' & n_pos_marker != 7 & k562 == 'cd19+']) + 
    geom_bar(aes(x=1, y=mean_pct, fill=n_pos_marker), stat='identity', width=1) + 
    facet_grid(t_type + day ~ car) + 
    scale_y_continuous(expand=c(0,0)) + 
    scale_x_discrete(expand=c(0,0)) +
    scale_fill_brewer('Number of Exh.\nMarkers Present', palette='YlOrRd') + 
    theme_minimal() +
    theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
    geom_vline(xintercept=0:7-0.5) +
    labs(x="CD4/CD8 and Day", y='CAR', title='Fraction of Cells (Donor 2)') +
    coord_polar("y", start=0)

Pie Charts, Mean

Remove day 33 because only donor 1.

marker_freq_avg <- marker_freq[day < 33, list(mean_pct= mean(mean_pct)), by=c('car', 'k562','n_pos_marker','t_type','day')]

ggplot(marker_freq_avg[k562 == 'cd19+']) + 
    geom_bar(aes(x=1, y=mean_pct, fill=n_pos_marker), stat='identity', width=1) + 
    facet_grid(t_type + day ~ car) + 
    scale_y_continuous(expand=c(0,0)) + 
    scale_x_discrete(expand=c(0,0)) +
    scale_fill_brewer('Number of Exh.\nMarkers Present', palette='YlOrRd') + 
    theme_minimal() +
    theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
    geom_vline(xintercept=0:7-0.5) +
    labs(x="CD4/CD8 and Day", y='CAR', title='Fraction of Cells (Donor Average)') +
    coord_polar("y", start=0)